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ABSTRACT 

In this paper we consider Basis Pursuit De-Noising (BPDN) 
problems in which the sparse original signal is drawn from 
a finite alphabet. To solve this problem we propose an iter- 
ative message passing algorithm, which capitalises not only 
on the sparsity but by means of a prior distribution also on 
the discrete nature of the original signal. In our numerical 
experiments we test this algorithm in combination with a 
Rademacher measurement matrix and a measurement matrix 
derived from the random demodulator, which enables com- 
pressive sampling of analogue signals. Our results show in 
both cases significant performance gains over a linear pro- 
gramming based approach to the considered BPDN problem. 
We also compare the proposed algorithm to a similar mes- 
sage passing based algorithm without prior knowledge and 
observe an even larger performance improvement. 

Index Terms — Compressive Sampling, Signal Recovery, 
Finite Alphabet, Message Passing 

1. INTRODUCTION 

When the information bearing part of a signal lies only in a 
small sub-space of the entire signal space, a uniform full rate 
sampling approach, i.e., sampling at Nyquist rate, is ineffi- 
cient. In 1 1 1 and [2|, Donoho and Candes et al. addressed 
this observation and introduced Compressive Sampling (CS). 
The key contribution of these papers was to show how random 
matrices and l\ -minimization can be applied to achieve opti- 
mal recovery of a sparse signal from a very limited number of 
measurements. 

While initial reconstruction algorithms were based on 
convex optimisation and linear programming, recent ad- 
vances in compressive sampling have led to the development 
of various other algorithms which solve the reconstruction 
task at a lower computational complexity. Inspired by the 
success of message passing algorithms as used for the de- 
coding operation of some channel codes, the authors of (3) 
solve the reconstruction of the original signal by means of 
belief propagation on a sparse graph. More recently, a simple 
iterative soft thresholding reconstruction algorithm, called 
Approximate Message Passing (AMP), has been proposed in 



0. The basis of this algorithm is also belief propagation, 
albeit on a fully connected graph, and it exhibits a virtually 
equivalent sparsity-undersampling trade-off to that of linear 
programming based reconstruction algorithms. 

In this paper, we address the problem of reconstructing a 
sparse finite alphabet signal from a limited number of noisy 
measurements which are obtained by compressive sampling. 
This problem appears in many areas such as spectrum sens- 
ing, symbol detection in digital communications, and multi- 
user detection, cf. ||5]|6). Many existing reconstruction algo- 
rithms for compressive sampling exploit the knowledge of the 
sparsity level of the original signal. Building upon the AMP 
framework, we propose here a novel AMP-based algorithm 
which capitalizes not only on the sparsity but also on a prior 
distribution, which manifests the finite alphabet property of 
the original signal. 

In the derivation of this algorithm we assume that a 
time discrete sparse signal is applied directly to a mea- 
surement matrix whose entries are randomly sampled from 
| — ^-7=, + -i= j, where R is the number of measurements. In 
practice however, we are typically confronted with various 
types of measurement matrices. To this end we consider here 
as an application also the random demodulator [7], which 
allows to sample time-continuous analogue signals at sub- 
Nyquist rates and resorts to techniques from compressive 
sampling to reconstruct the signal. In both cases our numeri- 
cal experiments indicate that the proposed algorithm offers an 
excellent performance in reconstructing sparse signals from 
noisy undersampled observations. When we compare the 
new algorithm with a standard linear programming based re- 
construction algorithm and the AMP algorithm for the BPDN 
problem from (8), we observe a significant performance im- 
provement. Moreover, we note that the computational cost of 
the algorithm proposed in this paper is comparable to existing 
state of the art algorithms. 

The remainder of the paper is organized as follows. In 
Section [2] we outline the considered compressive sampling 
problem. Next, we introduce our proposed approximate mes- 
sage passing reconstruction algorithm with prior knowledge 
in Section [3] In Section 14.11 we compare the proposed al- 
gorithm at first with the linear programming based SPGL1 



algorithm and the AMP algorithm for the BPDN problem in 
combination with a Rademacher measurement matrix. Then 
we briefly recap in Section 14.21 on the random demodulator 
and also present results for this scenario. Finally, Section [5] 
concludes the paper. 

2. PROBLEM OUTLINE 



argues that an estimate of the input distribution can be used 
to improve the recovery algorithm. Therefore, it is natural to 
employ a discrete prior distribution when the original signal 
is drawn from a finite alphabet, and we take 
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Let A = {ai, . . . , as} denote a finite set of S non-zero real 
numbers and define Ao := A U {0}. The VF-dimensional 
sparse column vector b shall have only K <C W non-zero 
entries, which are drawn uniformly at random from A. More- 
over, let * e R RxW be a matrix, for which R < W holds. 
Given the noisy observation 



*b 



(1) 



where n = [no, . . . , n r , . . . , n,R~i] T and n r is i.i.d. ]\f (0, er 2 ), 
our aim is to recover b. Problems of this type are frequently 
considered in the CS literature. One way to reconstruct the 
sparse vector b is by solving the optimization problem 



arg mm 



subject to 



*b 



< 



(2) 



for an optimization constant 7', which is chosen depending 
on the noise variance. In general though the Iq -minimization 
used in (01 is NP-hard. This has led to various approximate al- 
gorithms, including those based on convex relaxation, i.e., on 
ii-minimization, like BPDN H) or LASSO (TO). Many of the 
approximate algorithms for CS can thus be applied directly to 
approximate the solution of (f2]) by solving the relaxed prob- 
lem 



arg mm 



subject to 



*b 



< 



7, 



(3) 



and ensuring by, for example, thresholding that b g A^ 
holds. However, in addition to the number of non-zero entries 
in the sparse VK-dimensional signal b, the knowledge that its 
entries lie in the set Ao can also be utilized directly to im- 
prove the recovery of b, as will be discussed in the following 
sections. 

3. AMP WITH DISCRETE PRIOR DISTRIBUTION 

It is well known that the solution of the l\ -minimization in 
(0 corresponds to a mode of the posterior distribution when 
a double-exponential prior distribution is used. Donoho et 
al. also use a double-exponential prior in the derivation of 
their AMP algorithm [4], where the contribution of the prior 
distribution in the message update rules of the belief propa- 
gation algorithm can be interpreted as a 'sparsity promoting' 
soft thresholding operation. However, their results are quite 
general and a similar approach can be applied when a differ- 
ent choice of prior distribution is more suitable. Indeed, OD 



as prior for each b w in b = [bo, ...,b w , ...,bw-i] , where 
7To = 1 — ^7, and 7r s = for 1 < s < S. Note that this 
prior is constructed under the assumption that the non-zero 
entries in b are drawn uniformly from A. Should some non- 
zero entries be more likely than others, the prior distribution 
can be easily modified to reflect this additional information. 

Consider a fully connected bipartite graph between R 
measurement nodes on one side and W variable nodes on the 
other side. The measurement nodes shall represent the entries 
in v and likewise the variable nodes the unknown entries 
in b. As in [4 1 we study belief propagation (BP) message 
updates between the measurement and variable nodes on this 
complete graph. We denote the set of the W variable nodes 
as [W] and use w, uj € {0, 1, W — 1} as indices for this 
set. Similarly, we apply r, p G {0, 1, R — 1} as indices 
for the set [R] of all R measurement nodes. At iteration t we 
denote the message passed from the variable node w to the 
measurement node r by VwKr {b w ) and the message passed 
on this edge in the opposite direction by Vr% w {b w ). In the 
i-th iteration the BP message updates are then given as 
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where dh^ w denotes that integration is over all variables ex- 
cept b w , and a 2 is the variance of the noise in the measure- 
ments. Note that the messages from variable to measurement 
nodes in (0 are proportional to the probability mass function 
on _4o- We initialize 
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(7) 



We denote the mean and the variance of the message in © by 



QuUr and 



respectively. Note that £i°lj. r and Tw\ 



are 



initialized to the prior mean and variance Vw € [W] , W € 
[R] . We will derive the algorithm here under the assumption 

that -0 r u) £ | — ^7=, +^7t| j- Later in our simulations how- 
ever we relax this assumption and work with a general with 
normalized columns. 

Consider the random vector h_ w = [bo, 61, ... , b w _i, 
b w +i, . . . , b\Y-i] T distributed according to the product mea- 



n 
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(b^), and the associated scalar random 
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Fig. 1: Detection error rate P (b w b w j versus the noise variance <r 2 when 
the measurement matrix in (TJ is a Rademacher matrix with unit column 
norm, or the random demodulator matrix. The number of iterations for both 
AMP algorithms is set to T = 50. 



variable 
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Denote the induced density of xi% w by g. By the central 
limit theorem, for large W, g can be approximated by a Gaus- 



sian density with mean v r — 



UJ^W 



(*) 



and variance 



~k H/uj^w T u-+r- Thereby, if we write © in its equivalent 
form 
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we see that the factor to variable message update © can be 
approximated by a Gaussian integral. These observations lead 
to the simplified algorithm for message passing with a dis- 
crete prior distribution given in Algorithm [30. Note that we 
track the posterior probabilities in the log-domain, as numeri- 
cal simulations indicate that this results in a numerically more 
stable implementation of the simplified message-passing al- 
gorithm. 

4. NUMERICAL EXPERIMENTS 



a Rademacher matrix with unit column norm. To obtain the 
estimate b for the original b from the noisy observation v we 
apply the approximate message passing algorithm, as outlined 
in Algorithm [TJ For a comparison we also apply to the same 
problem the algorithms for the BPDN problem from ifTTlfTZll . 
named SPGL1-BPDN, and from [8), named AMP-BPDN. 
Unlike the proposed AMP algorithm the SPGL1-BPDN and 
the AMP-BPDN algorithm require the parameter 7 as input, 

which we choose as 7 = y/\W/R~\ - a-\fR- ^1 + 
In contrast to our proposed algorithm, which always returns 
an estimate for b in {-1, + 1} W , the SPGL1-BPDN and 
the AMP-BPDN algorithm return estimates in R w . For 
this reason one can threshold the b obtained from these two 
algorithms by a and take as the output sign ^5,^, where 

> a, and zero otherwise. Alternatively, one can search 

for the K entries in b with the largest magnitude and set them 
depending on their sign to ±1 and all other entries to zero. 

In FigureUJthe detection error rate, i.e., P (b w ^ b w ^j, is 

plotted for AMP-BPDN, SPGL1-BPDN, and the algorithm 
proposed in this paper versus the noise variance. Our sim- 
ulations indicate that the detection error rate of the SPGL1- 
BPDN and AMP-BPDN algorithm are strongly dependent 
on the chosen threshold a and that the decoding rule which 
simply chooses the K largest entries of b achieves the best 
performance for both of these algorithms (this version is plot- 
ted in Figure [TJ. However, even in this case the detection 
error probability performance of the SPGL1-BPDN and the 
AMP-BPDN algorithm is approximately 6dB respectively 
20dB worse than that of the proposed AMP algorithm with 
discrete prior. 

4.2. Random Demodulator 

In this section we consider discrete multi-tone signals, which 
occur, for example, in orthogonal frequency division multi- 
plex systems, in combination with the random demodulator. 
For this class, it is shown in |7| how the operation of the ran- 
dom demodulator, which is depicted in Figure [2] on the ana- 
logue signal can be described equivalently by a time-discrete 
representation. 

Let F £ £Wxw k e a (jiscj-gtg Fourier transform matrix 
and b £ { — 1,0, +1}^ a X-sparse vector. A time discrete 
representation of the analogue multi-tone signal is then given 



4.1. Rademacher Measurement Matrix 

For the numerical experiments presented here we assume that 
A = {-1,0, +1} and fix W = 512, R = 205, K = 20. 
The entries in the measurement matrix ^ shall be drawn uni- 
formly at random from | — l/\/R, +1/a/r|, which makes 



analogue 
input 
signal 




sub-Nyquist sampling rate 

compressed 
discrete 
signal 



chipping sequence 



1 Throughout the algorithm, "dot" is a placeholder for any r £ [R] or for 
0, in which case the summation in Step (3) is over all p £ [R]. 



Fig. 2: Structure of the random demodulator as discussed in [ 7 ] 



by x = Fb. The multiplication operation of the signal x with 
the chipping sequence and the integrate and dump operation 
of the random demodulator are described by the WxW diag- 
onal matrix D and by H e {0, 1} R xW , respectively. H shall 
have Wj B! consecutive ones in the ?'-th row starting from col- 
umn rW/R' + 1, where r = 0, . . . , R'— 1. As proposed in 
we allow H to have fractional elements in some of its columns 
when R does not divide W. To summarize, the sparse signal 
b is observed by the random demodulator through the random 



demodulator matrix 1 3>' = HDF e 



-_R.' x W 



. The sign of the 



non-zero entries of D on the main diagonal is chosen inde- 
pendently at random. However, we choose the magnitude of 
these entries here such that the columns of \P' have unit norm. 

As in Section |4~T1 we aim here to recover b from the noisy 
samples v given in (|T), where however the real valued mea- 
surement matrix 

"Re{*'}" 
Im{*'} 



(9) 



depends now on the random demodulator matrix \f. For our 
numerical experiments with the measurement matrix given by 
©, we set W = 512, R = 204, R' = 102, and K = 20. For 
this scenario the detection error rate, i.e., P (b w ^ b w ^j , is 
also plotted in Figure Q] versus the noise variance. The per- 
formance of all three algorithms is very similar to the case in 
Section l4~T1 with a Rademacher measurement matrix. The per- 
formance gap between the SPGL1-BPDN, the AMP-BPDN, 
and the proposed AMP algorithm with discrete prior is here 
also approximately 6dB and 20dB wide. 

5. CONCLUSIONS 



Algorithm 1 AMP with discrete prior 



• Input: Measurement matrix "J/, observation v, noise variance a 2 , 
alphabet Ao > prior probabilites 7To, 7Ti, . . . , tts, number of iterations 
T 

• Output: signal estimate b 

1. Initialize: t = 1, and l^l^ r (a s ) = tt s , for w £ [W], r £ [R], 

s e [S]. 

2. Calculate the mean and variance of the variable-to-factor messages: 

S 



s = 



t' w L+. = a s exp (a s ) - If. 

s=0 

3. Approximate the mean and variance of EIp^. ^p\w (b-w) 
= ^^P,™ \ V P ~ X} 



'/ 



1 

R 



A*) 



UJ^W 



4. Incorporate prior and normalize: 

f£U. (a s ) = log?r s - - 



(*) _ 



2 V 
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AD 



(a s ) 



log I ^2 cxpl^K. (a s >) ] +lwh. («s) 



Set t 4— t + 1, and repeat (2)-(4) until stopping criterion holds. 



5. Return b = \bo, ■ ■ ■ , bw—i , where 



b w = arg max l^L,* (a) 



In this paper we have developed a novel reconstruction algo- 
rithm for compressive sampling problems, which applies the 
prior knowledge that the entries of the original signal vec- 
tor belong to a finite alphabet. Our simulation results show 
for this algorithm significant performance gains over existing 
reconstruction algorithms for the BPDN problem. In future 
work we hope to further simplify the message passing update 
equations of the proposed algorithm and thus further reduce 
its complexity. 
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